Toward a Universal Formulation of the Halo Mass Function 
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We compute the dark matter halo mass function using the excursion set formalism for a diffusive 
barrier with linearly drifting average which captures the main features of the ellipsoidal collapse 
model. We evaluate the non-Markovian corrections due to the sharp filtering of the linear density 
field in real space with a path-integral method. We find an unprecedented agreement with N-body 
simulation data with deviations < 5% over the range of masses probed by the simulations. This 
indicates that the Excursion Set in combination with a realistic modelling of the collapse threshold 
can provide a robust estimation of the halo mass function. 
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A large body of evidence suggests that dark matter 
(DM) plays a crucial role in the formation, evolution and 
spatial distribution of cosmic structures \^tM ■ Central to 
the DM paradigm is the idea that initial density fluctua- 
tions grow under gravitational instability eventually col- 
lapsing into virialized objects, the halos. It is inside these 
gravitationally bounded structures that cooling baryonic 
gas falls in to form the stars and galaxies we observe 
today. Consequently, the study of the halo mass distri- 
bution is of primary importance in cosmology. In the 
Press-Schechter approach Q , the number of halos in the 
mass range [M, M + dM] can be written as 



dn _ p dlogg ^ 

dM ~ ■'^^'Jr^ dlogM 



(1) 



where p is the background matter density and a{M) is 
the root-mean-square fluctuation of the linear dark mat- 
ter density fleld smoothed on a scale R{M) (containing 
a mass M), with 



dk k^P{k)W^[k,R{M)], (2) 



where P{k) is the linear DM power spectrum and 
W{k, R) is the Fourier transform of the smoothing (fil- 
ter) function in real space. In Eq. ([l}, the function 
/((t) — 2a^F{a^), known as 'multiplicity function', en- 
codes the effects of the gravitational processes responsi- 
ble for the formation of halos through its dependence on 
F{S) = dF/dS, with F{S) being the fraction of mass 
elements in halos of mass > M{S). Hereafter, we will 
refer to /(cr) simply as the halo mass function. 

The collapse of halos is a highly nonlinear gravitational 
process that has been primarily investigated using nu- 
merical N-body simulations. Over the past few years 
several numerical studies have measured /(cr) at few per- 
cent uncertainty level for various cosmologies and using 
different halo detection algorithms (see e.g. [§-11]). On 
the other hand, we still lack an accurate theoretical esti- 
mation of the halo mass function. Following the seminal 



work by Press and Schechter the excursion set the- 
ory 10] has provided us with a consistent mathematical 
framework for computing /(u) from the statistical prop- 
erties of the initial DM density field (for a review see [11|). 
Nevertheless, an analytical derivation of /(cr) can be ob- 
tained only for a top-hat filter in Fourier space (sharp-k 
filter). Although Monte-Carlo simulations can be used 
in the case of generic filters (see e.g. 13, 121), most of 
the work in the literature has focused on the modeling of 
the halo collapse conditions and the comparison with N- 
body simulations, while assuming the sharp-fc filter (see. 



e.g., [I3|-|l6j). However, such a smoothing function does 



not correspond to any realistic halo mass definition. The 
issue has been recently addressed by Maggiore and Ri- 
otto who made a major contribution by introducing 
a path-integral method that extends the analytical com- 
putation to generic filters. 

In this Letter we present the first thorough comparison 
against N-body simulation data of the excursion set mass 
function with top-hat filter in real space for a stochastic 
barrier model which encapsulates the main characteris- 
tics of the ellipsoidal collapse of dark matter. A detailed 
derivation of these results is given in a companion paper 

Let us consider the DM density contrast, (5(x), 
smoothed on the scale R, 



<5(x,i?) = J d^yW{\^-y\,R)Siy), 



(3) 



where M^(x, K) is the smoothing function in real space. 
Bond et al. [lOj have shown that at any given point in 
space, 5(x, R) performs a random walk as a function of 
the variance of the smoothed linear density field S{R). 
The formation of halos of mass M corresponds to tra- 
jectories S{S) crossing for the first time a barrier B at 
S{M), i.e., 6{S) = B, where the value of B depends 
on the assumed gravitational collapse criterion. In the 
case of the spherical collapse model [l^ B — 6c-, that is 
the linearly extrapolated density of a top-hat spherical 



2 



perturbation at the time of collapse. Then, the evalua- 
tion of f{a) is reduced to computing the rate at which 
the random walks hit the barrier for the first time, i.e., 
T{S) = dF/dS. 

The nature of the random walk depends on the filter- 
ing procedure, which specifies the relation between the 
smoothing scale R and the halo mass definition M. For 
a sharp-fc filter, W{k,R) = 9{1/R— fc), and Gaussian 
initial conditions, S{S) performs a Markov random walk 
described by the Langevin equation: 

^=^siS), (4) 

with noise 775(5*) such that (775(6')) = and 
{risiS)r]s{S')) = doiS - S'), where So is the Dirac- 
function (for the full derivation, see, e.g., [ll|, 17|). As 
first shown in , the probability distribution of the tra- 
jectories satisfies a simple Fokker-Planck equation with 
absorbing boundary at 6{S) = Sc- The resulting first- 
crossing distribution gives the Press-Schechter formula 
with the correct normalization factor (the so called 
'extended Press-Schechter'). 

However, the spherical collapse model is a simplistic 
approximation of the nonlinear evolution of matter den- 
sity fiuctuations. As shown in po| . initial Gaussian per- 
turbations are highly nonspherical. Hence, the collapse 
of a homogeneous ellipsoid (see, e.g., [21|) should pro- 
vide a far better description. In such a model the critical 
density threshold depends on the eigenvalues of the de- 
formation tensor, which are random variables with prob- 
ability distributions _that d epend o n the statistics of the 



linear density field [14 



Because of this, the 
barrier behaves as a stochastic variable itself, performing 
a random walk whose properties depend on the speci- 
ficities of the collapse model considered. For example, 
Sheth et al. [3 showed that the average of the barrier is 
{B{S)) = Sc[l+l3iS/6l)'^], with 13 = 0.47 and 7 = 0.615. 

The recent analysis of halos in N-body simulations has 
confirmed the stochastic barrier hypothesis 27[. Mag- 



giore and Riotto [28[ have modeled these features assum- 
ing a stochastic barrier with average {B{S)) = Sc and 
variance {{B — {B{S)))^) = S Db, where Db is a con- 
stant diffusion coefficient. Here, we improve their barrier 
model by assuming a Gaussian diffusion with linearly 
drifting average {B{S)) — Sc + [l3[ which approxi- 
mates the ellipsoidal collapse prediction [ij] . Recently, a 
general analysis of nondiffusive moving barriers has been 
presented in (29j . However, this work has mainly focused 
on the mass function in the presence of Non-Gaussian 
initial conditions rather than the comparison with Gaus- 
sian N-body simulations. The Langevin equation for this 
barrier model reads as 



— =/3 + ,7b(5), 



(5) 



where the noise 773(5) is characterized by (773(6')) — 
and (77b(6')77b(6')) = DbSd{S - S'). Without loss of 



generality we can assume that t]b{S) and r]s{S) are un- 
correlated. It is convenient to introduce Y = B — S and 
rewrite Eqs. ^ and ([5]) as a single Langevin equation: 



dY 



(6) 



with white noise r;(6) = 775(6")-)- 7^3(6') such that (77(6)) — 
and (77(6')?7(6")) = (1 + Db)S{S ~ S'). The Fokker- 
Planck equation associated with Eq. ^ and describing 
the probability Ho(yo, S) reads as 



as 



dY 



1 + Db d^Ho 
2 dY^ 



(7) 



where we indicate with the "0" underscore the fact that 
Hg is associated to a Markov process. 

The system starts at {(5(0) = 0,-6(0) — (5c}; hence, 
we solve Eq. ([7]) with initial condition Yq = Sc and 
impose the absorbing boundary condition at Y = 0, 
i.e., Ho(0,6) = 0. For a concise notation we omit the 
dependence on Yq and simply refer to Hq {Y, S) . By 
rescaling the variable Y Y = Y j \J\ -I- Db^ a fac- 
torizable solution can be found in the form Ho(y, 6") = 
J7(y,6)exp[c(f - c6'/2)], where c = /3/Vl + -Ds and 
C/(F, 6) satisfies a diffusion equation. Using the above 
initial condition, the latter can be solved with the image 
method [30| or by Fourier transform. Thus, we obtain 



no(F,6) 



^2-kS(\ + Db) 



g 2S(1 + Db) 



g 2S(1 + Db) 



(8) 

In general the Fokker-Planck equation for random walks 
with nonlinear biased diffusion and absorbing boundary 
condition does not have an exact analytic solution. This 
is why we have assumed the linearly drifting average 
barrier rather than the prediction of the ellipsoidal col- 
lapse model IJ]. As we will see later, having an ex- 



act analytical solution greatly simplify the evaluation 
of the corrections due to the smoothing function. We 
should remark that the above solution is defined only for 
Y > Since the number of trajectories is conserved, 
then the first-crossing distribution is obtained by deriv- 
ing \^T^{S')dS' = 1 - J^no{Y,S)dY from which we 
finally obtain the Markovian mass function 



/o(fT) = 



(T^/l + Db 



2<t2(1 + d_b) 



(9) 



for /? = and Db = this coincides with the stan- 
dard Markovian solution that gives the extended Press- 
Schechter formula, while for = we recover the so- 
lution for the nondiffusive linearly drifting barrier 11|. 

As mentioned earlier, a crucial point of this deriva- 
tion is the assumption of the sharp-fc filter. In numer- 
ical N-body simulations the mass definition depends on 
the halo detection algorithm. For instance, the spherical 
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overdensity (SOD) halo finder detects halos as groups 
of particles in a spherical regions of radius i?A contain- 
ing a density pA = with A an overdensity param- 
eter usually fixed to A = 200. Thus, the halo mass is 
M = 4/37ri?^pA, which is equivalent to having a sharp- 
X filter, or W{k,R) = (3/fci?)[sin(fci?) - (kR) cos{kR)]. 
However, in this case the stochastic evolution of the sys- 
tem is no longer Markovian. Hence, in order to con- 
sistently compare the excursion set mass function with 
SOD estimates of /(cr) it is necessary to account for the 
correlations induced by W(k, R). 

Maggiore and Riotto [17| have shown that these cor- 
relations can be treated as perturbations about the 
"zero" -order Markovian solution. More specifically, the 
noise variable r]{S) acquires a perturbative correction, 
{r,{S)T]{S')) = {l+DB)SDiS-S')+AiS,S'), which in 
the case of the sharp-x filter can be approximated by 
A{S,S') « kS{S' ~ S)/S'. For the concordance A Cold 
DM model we find k « 0.47. Using the path-integral 
technique described in 17i], we compute the corrections 
to Ho(y, S) to first order in k. These consist of a "mem- 
ory" term, 

U^ = -dY / dS'AiS',S)uliYo,0,S')nl{0,Y,S-S'), 
Jo 

(10) 

and a "memory-of-memory" term 

H™-" = I^dS' Js„ dS"A{S',S")nl{YQ,0,S')x 

xH^(0,0,5"-5')H^(0,y,5-5'), (11) 

where H^(yo,0,S'), H^(0,y,5) and H^(0,0,5) in 
Eqs. pU)) and ((TT|) are given by the finite time corrections 
of the Markovian solution near the barrier (see [l3|)- We 
find 



where a 
lytically, we find 



n^(Fo,0,5) 




a(Y-o + /3S)2 


(12) 




e 26' 


HiJ(0,y,5) 


aY 




(13) 




-e 2s , 

TT 


H^(0,0,5) 






(14) 


^ 1/(1 + ^B. 


. Eq. (nni 


) can be computed 


ana- 



li^ = -kaYody 




yga^*(y-yo-/3|)Erfc ,, 

25 

(15) 

where k = n/ (I+Db)- Since Equation (1151) is linear in Y, 
the integration of T^{S) = -d/dS H™ dY vanishes. 
Thus, the memory term does not contribute to the mass 
function independently of the barrier behavior (in agree- 
ment with [III). The double integral in the memory-of- 
memory term cannot be computed analytically, in such a 
case we expand the integrands in powers of /? (given that 
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FIG. 1: Contributions to the halo mass function /tot (solid 
line) for /? — 0.2 and Db ~ 0.6. The different curves cor- 
respond to the Markovian mass function fo (dotted line), 
f^-^to (short-dashed line), /J^^-^t, (long-dashed line), /J^^-^,, 
(dot-short dashed line), Z™^™) (dot-long dashed line). 



from the ellipsoidal collapse we expect (3 < 1). By com- 
puting Tr~""{S) = -d/dSj^ n™-'"(iy and expressing 
the results directly in terms of /(cr), we find the non- 
Markovian correction to zero order in (3 (i.e. (3 — 0) to 
be 



fm — m { „\ 



e 2<7 



irfo^ 

2 r'2a2 



(16) 

where r(0, z) is the incomplete Gamma function. Not 
surprisingly this expression coincides with the memory- 
of-memory term in [ITj . The first order correction in /3 
is given by 



/!V=o(^)+^Erfc 



(17) 



and the second order reads 



{a)=P'a5,k\a Erfc i^J- 



a5l 



3a(52 



+T-fr 0,-^ 



4 fT2 



2a2 



(18) 



For /3/(l + Db) < 1, corrections 0{> P^) are negligible 
(see, e.g.. Fig. [T|), hence, Eqs. ® and ([Tn])-(ITH]) give the 
relevant contributions to the mass function. 
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FIG. 2: (Upper panel) Halo mass function a.t z — given by 
the Tinker et al. fitting formula for A = 200 (solid blue line), 
diffusing drifting barrier with j3 = 0.057 and Dj, — 0.294 (red 
dashed line) and Maggiore & Riotto [2^ with Db ~ 0.235 
(green dotted line). Data points are from [^. (Lower panel) 
Relative difference with respect to the Tinker et al. fitting 
formula. The thin black solid lines indicates 5% deviations. 



rors ajs = 0.001 and = 0.001 respectively. In Fig. [2] 
(upper panel) we plot the corresponding mass function 
(red dash line) against the simulation data together with 
the four-parameter fitting formula by Tinker et al. for 
A = 200 (solid blue line). For comparison we also plot 
the diffusive barrier case by Maggiore and Riotto (isj 
which best fit the data with Db — 0.235 (green dotted 
line). In Fig. [2] (lower panel) we plot the relative differ- 
ences with respect to the Tinker et al. formula. We may 
notice the remarkable agreement of the diffusive drifting 
barrier with the data. Deviations with respect to Tin- 
ker et al. (2008) are < 5% level over the range of masses 
probed by the simulations. This is quite impressive given 
the fact that our model depends only on two physically 
motivated parameters. 

In the upcoming years a variety of astrophysical ob- 
servations will directly probe dn/dM. The halo mass 
function we have derived here can provide the base for 
a through cosmological model comparison. In a com- 
panion paper we will describe in detail the derivation of 
these results, as well as extensive discussion on the red- 
shift evolution of the mass function and halo bias. 

We are especially thankful to J. Tinker for kindly pro- 
viding us with the mass function data. It is a pleasure 
to thank J.-M. Alimi, Y. Rascra, T. Riotto and R. Sheth 
for useful discussions. I. Achitouv is supported by the 
"Ministere de I'Education Nationale, de la Recherche et 
de la Technologie" (MENRT). 



In principle the values of /3 and Dg as well as their 
redshift and cosmology dependence can be predicted in 
a given halo collapse model by computing the average 
and variance of the probability distribution of the col- 
lapse density threshold. However, this requires a dedi- 
cated study which should also include environmental ef- 
fects that have been shown to play an important role in 
determining the properties of the halo mass distribution 
[2^ . This goes beyond the scope of this Letter. 

Here, we take a different approach. /3 and Db are 
physical motivated model parameters which we can cali- 
brate against N-body simulation data, and test whether 
the mass function derived above provides an acceptable 
description of the data. To this purpose we use the mea- 
surements of the halo mass function obtained by Tin- 
ker et al. [i] using SOD(200) on a set of WMAP-1 yr 
and WMAP-3 yr cosmological N-body simulations. For 
these cosmological models the spherical collapse predicts 
5c = 1.673 at z = (for a detailed calculation see [l]). 
Using such a value, we run a likelihood Markov chain 
Monte Carlo analysis to confront the mass function previ- 
ously computed against the data at z = in the prior pa- 
rameter space log/3 = [—4,0] and logDg = [—3,0]. We 
find the best fit values to be /3 = 0.057 and Db = 0.294. 
The data strongly constrain these parameters, with er- 
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